Author

Tim Kerr

Code
### Packages
library(haven)
library(rstan)
library(loo)
library(StanHeaders)
library(tidyverse)
library(reshape2)
library(ggplot2)
library(ggExtra)
library(ggforce)
library(GGally)
library(patchwork)
library(gganimate)
library(MetBrewer)
library(ggh4x)
library(transformr)
library(data.table)
library(coda)
library(reactable)
library(kableExtra)
# library(reactablefmtr)
library(car)
library(mvtnorm)
library(ggdist)
library(tidybayes)
# library(rstatix)
library(easystats)
library(tidymodels)
library(broom)
library(pixiedust)
library(broom.mixed)

library(brms)
library(cmdstanr)
library(posterior)
library(bayesplot)
Code
cp.cs <- met.brewer("Hiroshige")[c(1,10,2,8)]
cp.cs <- append(cp.cs, met.brewer("Cassatt2")[c(3,8)])

theme_tim <- function(base_size = 14) {
  theme_modern(base_size = base_size) %+replace%
    theme(
      # L'ensemble de la figure
      plot.title = element_text(size = rel(1), face = "bold", margin = margin(0,0,5,0), hjust = 0),
      # Zone où se situe le graphique
      panel.grid.minor = element_blank(),
      panel.grid.major = element_blank(),
      panel.border = element_blank(),
      # Les axes
      axis.title = element_text(size = rel(0.85), face = "bold"),
      axis.text = element_text(size = rel(0.70), face = "bold"),
      axis.line = element_line(color = "black",),
      # La légende
      legend.title = element_text(size = rel(0.85), face = "bold"),
      legend.text = element_text(size = rel(0.70), face = "bold"),
      legend.key = element_rect(fill = "transparent", colour = NA),
      legend.key.size = unit(1.5, "lines"),
      legend.background = element_rect(fill = "transparent", colour = NA),
      # Les étiquettes dans le cas d'un facetting
      strip.background = element_rect(fill = "#17252D", color = "#17252D"),
      strip.text = element_text(size = rel(0.85), face = "bold", color = "white", margin = margin(5,0,5,0))
    )
} 
# Changing the default theme
theme_set(theme_tim(base_size = 10))
Code
### cmdStanr stuff

check_cmdstan_toolchain()

# install_cmdstan(cores = 4, overwrite = TRUE)

cmdstan_path()
[1] "/Users/timkerr/.cmdstan/cmdstan-2.33.1"
Code
cmdstan_version()
[1] "2.33.1"
Code
options(
  mc.cores = 4,
  brms.threads = 4,
  brms.backend = "cmdstanr",
  brms.file_refit = "on_change"
)
Code
merge_prior_post <- function(gen,rec,...){
df1 <- gen %>%
  gather_draws(...) %>% ungroup() %>% mutate(model = "ground truth")

df2 <- rec %>%
  gather_draws(...) %>% ungroup() %>% mutate(model = "posterior")

df <- bind_rows(df1,df2) %>% mutate(model = as.factor(model))

return(df)
}

Model spec

Data

Assuming a simple Bernoulli logit task / model, with individual/subject level theta determining outcomes on the task. (In reality I’m using a reinforcement learning task, but this is a simple approximation to start with).

\[ y \sim \text{BernoulliLogit}(\theta)\] I’m assuming/testing that a psychometric measure, like the GAD-7, has an impact on this task performance. I think the impact should be bi-directional rather than additive (i.e. low GAD scores improve task score, high GAD scores mean worse task performance)

GAD is distributed roughly exponentially at population level, most have scores between 0 and 5, then a long tail up to the maximum of 21.

GAD GAD

To facilitate the bi-directionality of effect, standardising around a mean of zero seems an OK way to start with this (would ideally estimate the fuzzy point at which the direction of effect would change).

\[\psi = \text{Exponential}(1)\]

\[z_{\psi} = \frac{\psi - \bar\psi}{\sigma_{\psi}} \]

Parameters

The aim then is to estimate the group level impact of GAD-7 ( \(\psi\) ) on the task, via the parameter \(m\), a scalar/weight between 0 and 1.

\[ p(\theta) = \text{Normal}(\Theta_{pop} + m\psi,\sigma_{pop}) \]

Data generator

I used Stan to generate data in one script, generating \(y\) from a Bernoulli rng, and \(\psi\) from an exponential rng. I inputted priors at set values to generate this data from these prior values.

\[\begin{aligned} &\large \bf Priors \\ \Theta^{Group/pop} &\sim \text{Normal}(0.3,0.1) \\ \sigma^{Group/pop} &\sim \text{Normal}(0.1,0.1) \\ m &\sim \text{Normal}(0.3,0.01) \\ \psi &\sim \text{Exponential}(1) \end{aligned}\]

Running the model

Recovery

Then inputted the generated data into the same model, and compared ground truth to posterior/recovered parameters.

Priors were uninformative:

\[\begin{aligned} &\large \bf Priors \\ \Theta^{Group/pop} &\sim \text{Normal}(0,1) \\ \sigma^{Group/pop} &\sim \text{Normal}(0,1) \\ m &\sim \text{Normal}(0,1) \\ \end{aligned}\]

Whilst it recovers the group level theta mean and sd parameters well, it sverely underestimates the \(m\) parameter, the one of interest.

I wonder if this is due to the exponential distribution of GAD, given most people’s scores will sit in the bulk of the distribution and not have much of an effect. Perhaps this score needs transforming on the log scale or suchlike to make this work…

Code
file <- "~/FLARe/FLARe_Modelling/Stan_experimental/Embedding_Generator.stan"

mod <- cmdstan_model(file)

mod$check_syntax(pedantic = TRUE)

t_mult <- 1

datalist <- list(N = 100,
                 T = 10*t_mult,
   gp_theta_mu_pr = 0.3,
   gp_theta_sigma_pr = 0.1, 
   gp_sigma_mu_pr = 0.1,
   gp_sigma_sigma_pr = 0.1,
  gp_m_mu_pr = 0.3,
  gp_m_sigma_pr = 0.01,
  psi_pr = 1
)

fit_gen <- mod$sample(
  data = datalist,
  seed = 123,
  chains = 4,
  parallel_chains = 4,
  refresh = 1000,
  iter_warmup = 250,
  iter_sampling = 1000,
  # fixed_param = TRUE
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 2 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 3 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 4 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 1 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 2 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 3 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 4 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 1 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 1 finished in 1.3 seconds.
Chain 2 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 4 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 2 finished in 2.0 seconds.
Chain 4 finished in 2.0 seconds.
Chain 3 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 3 finished in 3.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.1 seconds.
Total execution time: 3.2 seconds.
Code
theta <- fit_gen %>% spread_draws(subj_theta[n]) %>% sample_draws(ndraws = 1, seed = 789) %>% pull(subj_theta)
psi <- fit_gen %>% spread_draws(psi_pred[n]) %>% sample_draws(ndraws = 1, seed = 789) %>% pull(psi_pred)
y <- fit_gen %>% spread_draws(y_pred[n,t]) %>% sample_draws(ndraws = 1, seed = 789) %>% pull(y_pred) %>% matrix(nrow = 10*t_mult, ncol =100) %>% t()

file <- "~/FLARe/FLARe_Modelling/Stan_experimental/Embedding_Recovery.stan"

mod <- cmdstan_model(file)

mod$check_syntax(pedantic = TRUE)

datalist <- list(N = 100,
                 T = 10*t_mult,
                 y = y,
   gp_theta_mu_pr = 0,
   gp_theta_sigma_pr = 1, 
   gp_sigma_mu_pr = 0,
   gp_sigma_sigma_pr = 1,
  gp_m_mu_pr = 0,
  gp_m_sigma_pr = 1,
  psi_y = psi
)

fit_rec <- mod$sample(
  data = datalist,
  seed = 123,
  chains = 4,
  parallel_chains = 4,
  refresh = 1000,
  iter_warmup = 250,
  iter_sampling = 1000,
  
  # fixed_param = TRUE
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 2 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 3 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 4 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 1 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 2 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 3 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 4 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 3 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 3 finished in 3.4 seconds.
Chain 4 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 4 finished in 3.4 seconds.
Chain 2 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 2 finished in 3.8 seconds.
Chain 1 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 1 finished in 4.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 3.7 seconds.
Total execution time: 4.2 seconds.
Code
summ <- fit_gen$summary()

summ[1:4,] %>% kable(digits = 3, caption = "Ground Truth")
Ground Truth
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -24.550 -26.295 43.273 46.404 -94.294 49.958 1.026 97.118 168.323
gp_theta 0.298 0.299 0.099 0.101 0.138 0.458 1.003 809.861 1548.767
gp_sigma 0.165 0.157 0.067 0.069 0.073 0.289 1.024 100.588 189.324
gp_m 0.300 0.300 0.010 0.010 0.284 0.317 1.002 8013.911 2918.352
Code
summ2 <- fit_rec$summary()

summ2[1:4,] %>% kable(digits = 3, caption = "Recovery")
Recovery
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -620.137 -627.061 33.713 25.027 -659.058 -546.834 1.047 69.719 42.226
gp_theta 0.156 0.157 0.078 0.078 0.025 0.286 1.012 353.531 245.667
gp_sigma 0.363 0.366 0.118 0.114 0.148 0.553 1.043 74.198 43.094
gp_m 0.095 0.089 0.060 0.061 0.012 0.205 1.004 772.950 909.604
Code
merged_df <- merge_prior_post(fit_gen,fit_rec,gp_theta, gp_sigma,gp_m)

merged_df %>% 
  # filter(model == "prior") %>% 
 ggplot(aes(y = .variable, x = .value)) + 
  stat_dotsinterval(aes(group = model, 
                        slab_colour = model, 
                        slab_fill = model,
                        interval_colour = model,
                        point_colour = model)) +
  facet_wrap(~.variable, scales = "free_x") +
  ggtitle("PPC of group level parameters")

Code
#ppc
fit_rec %>% spread_draws(subj_theta[n]) %>% ggplot(aes(x = subj_theta)) + stat_dots(quantiles = 100, fill = "red", alpha = 0.5) + stat_dots(data = as.data.frame(theta), aes(x = theta), quantiles = 100, fill = "blue", alpha = 0.5, inherit.aes = TRUE) + ggtitle("PPC of subject level theta parameters, Red = generated/estimated, Blue = real/ground truth")

Code
fit_rec %>% spread_draws(y_pred[n,t]) %>% sample_draws(ndraws = 1, seed = 789) %>% ggplot(aes(x = y_pred)) + stat_dots(quantiles = 100, fill = "red", alpha = 0.5) + stat_dots(data = as.data.frame(as.vector(y)), aes(x = `as.vector(y)`+0.1), quantiles = 100, fill = "blue", alpha = 0.5, inherit.aes = TRUE) + ggtitle("PPC Y, Red = generated/estimated, Blue = real/ground truth") + ggtitle("PPC of Bernoulli data, Red = generated/estimated, Blue = real/ground truth (shifted right to see easier)")

Changing the priors

Running the same but with more precise/informative priors

\[\begin{aligned} &\large \bf Priors \\ \Theta^{Group/pop} &\sim \text{Normal}(0,1) \\ \sigma^{Group/pop} &\sim \text{Normal}(0,1) \\ m &\sim \text{Normal}(0.3,0.1) \\ \end{aligned}\]

Code
file <- "~/FLARe/FLARe_Modelling/Stan_experimental/Embedding_Recovery.stan"

mod <- cmdstan_model(file)

mod$check_syntax(pedantic = TRUE)

datalist <- list(N = 100,
                 T = 10*t_mult,
                 y = y,
   gp_theta_mu_pr = 0,
   gp_theta_sigma_pr = 1, 
   gp_sigma_mu_pr = 0,
   gp_sigma_sigma_pr = 1,
  gp_m_mu_pr = 0.3,
  gp_m_sigma_pr = 0.1,
  psi_y = psi
)

fit_rec <- mod$sample(
  data = datalist,
  seed = 123,
  chains = 4,
  parallel_chains = 4,
  refresh = 1000,
  iter_warmup = 250,
  iter_sampling = 1000,
  
  # fixed_param = TRUE
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 2 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 3 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 4 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 2 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 3 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 1 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 4 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 4 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 4 finished in 3.8 seconds.
Chain 2 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 3 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 2 finished in 4.3 seconds.
Chain 3 finished in 4.2 seconds.
Chain 1 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 1 finished in 5.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 4.4 seconds.
Total execution time: 5.3 seconds.
Code
summ <- fit_gen$summary()

summ[1:4,] %>% kable(digits = 3, caption = "Ground Truth")
Ground Truth
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -24.550 -26.295 43.273 46.404 -94.294 49.958 1.026 97.118 168.323
gp_theta 0.298 0.299 0.099 0.101 0.138 0.458 1.003 809.861 1548.767
gp_sigma 0.165 0.157 0.067 0.069 0.073 0.289 1.024 100.588 189.324
gp_m 0.300 0.300 0.010 0.010 0.284 0.317 1.002 8013.911 2918.352
Code
summ2 <- fit_rec$summary()

summ2[1:4,] %>% kable(digits = 3, caption = "Recovery")
Recovery
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -625.027 -631.354 29.510 24.954 -661.374 -565.642 1.016 138.622 117.134
gp_theta 0.162 0.163 0.074 0.073 0.040 0.284 1.001 1299.264 2234.286
gp_sigma 0.376 0.379 0.119 0.123 0.172 0.566 1.016 144.997 130.201
gp_m 0.161 0.160 0.062 0.062 0.058 0.263 1.007 1065.776 1066.513
Code
merged_df <- merge_prior_post(fit_gen,fit_rec,gp_theta, gp_sigma,gp_m)

merged_df %>% 
  # filter(model == "prior") %>% 
 ggplot(aes(y = .variable, x = .value)) + 
  stat_dotsinterval(aes(group = model, 
                        slab_colour = model, 
                        slab_fill = model,
                        interval_colour = model,
                        point_colour = model)) +
  facet_wrap(~.variable, scales = "free_x") +
  ggtitle("PPC of group level parameters")

Code
#ppc
fit_rec %>% spread_draws(subj_theta[n]) %>% ggplot(aes(x = subj_theta)) + stat_dots(quantiles = 100, fill = "red", alpha = 0.5) + stat_dots(data = as.data.frame(theta), aes(x = theta), quantiles = 100, fill = "blue", alpha = 0.5, inherit.aes = TRUE) + ggtitle("PPC of subject level theta parameters, Red = generated/estimated, Blue = real/ground truth")

Code
fit_rec %>% spread_draws(y_pred[n,t]) %>% sample_draws(ndraws = 1, seed = 789) %>% ggplot(aes(x = y_pred)) + stat_dots(quantiles = 100, fill = "red", alpha = 0.5) + stat_dots(data = as.data.frame(as.vector(y)), aes(x = `as.vector(y)`+0.1), quantiles = 100, fill = "blue", alpha = 0.5, inherit.aes = TRUE) + ggtitle("PPC Y, Red = generated/estimated, Blue = real/ground truth") + ggtitle("PPC of Bernoulli data, Red = generated/estimated, Blue = real/ground truth (shifted right to see easier)")

Increasing the magnitude of effect

Is the effect too subtle to be picked up? Increasing the value of \(m\) to 0.7…

Code
file <- "~/FLARe/FLARe_Modelling/Stan_experimental/Embedding_Generator.stan"

mod <- cmdstan_model(file)

mod$check_syntax(pedantic = TRUE)

t_mult <- 1

datalist <- list(N = 100,
                 T = 10*t_mult,
   gp_theta_mu_pr = 0.3,
   gp_theta_sigma_pr = 0.1, 
   gp_sigma_mu_pr = 0.1,
   gp_sigma_sigma_pr = 0.1,
  gp_m_mu_pr = 0.7,
  gp_m_sigma_pr = 0.01,
  psi_pr = 1
)

fit_gen <- mod$sample(
  data = datalist,
  seed = 123,
  chains = 4,
  parallel_chains = 4,
  refresh = 1000,
  iter_warmup = 250,
  iter_sampling = 1000,
  # fixed_param = TRUE
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 2 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 3 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 4 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 1 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 3 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 2 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 4 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 1 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 1 finished in 2.1 seconds.
Chain 3 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 3 finished in 2.3 seconds.
Chain 2 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 2 finished in 3.0 seconds.
Chain 4 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 4 finished in 3.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.7 seconds.
Total execution time: 3.6 seconds.
Code
theta <- fit_gen %>% spread_draws(subj_theta[n]) %>% sample_draws(ndraws = 1, seed = 789) %>% pull(subj_theta)
psi <- fit_gen %>% spread_draws(psi_pred[n]) %>% sample_draws(ndraws = 1, seed = 789) %>% pull(psi_pred)
y <- fit_gen %>% spread_draws(y_pred[n,t]) %>% sample_draws(ndraws = 1, seed = 789) %>% pull(y_pred) %>% matrix(nrow = 10*t_mult, ncol =100) %>% t()

file <- "~/FLARe/FLARe_Modelling/Stan_experimental/Embedding_Recovery.stan"

mod <- cmdstan_model(file)

mod$check_syntax(pedantic = TRUE)

datalist <- list(N = 100,
                 T = 10*t_mult,
                 y = y,
   gp_theta_mu_pr = 0,
   gp_theta_sigma_pr = 1, 
   gp_sigma_mu_pr = 0,
   gp_sigma_sigma_pr = 1,
  gp_m_mu_pr = 0,
  gp_m_sigma_pr = 1,
  psi_y = psi
)

fit_rec <- mod$sample(
  data = datalist,
  seed = 123,
  chains = 4,
  parallel_chains = 4,
  refresh = 1000,
  iter_warmup = 250,
  iter_sampling = 1000,
  
  # fixed_param = TRUE
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 2 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 3 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 4 Iteration:    1 / 1250 [  0%]  (Warmup) 
Chain 1 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 2 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 3 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 4 Iteration:  251 / 1250 [ 20%]  (Sampling) 
Chain 1 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 2 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 1 finished in 3.5 seconds.
Chain 2 finished in 3.5 seconds.
Chain 3 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 4 Iteration: 1250 / 1250 [100%]  (Sampling) 
Chain 3 finished in 3.5 seconds.
Chain 4 finished in 3.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 3.5 seconds.
Total execution time: 3.7 seconds.
Code
summ <- fit_gen$summary()

summ[1:4,] %>% kable(digits = 3, caption = "Ground Truth")
Ground Truth
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -20.677 -22.029 42.254 42.511 -87.460 52.369 1.039 53.852 78.233
gp_theta 0.302 0.300 0.096 0.093 0.142 0.464 1.003 1611.850 1280.224
gp_sigma 0.159 0.150 0.063 0.063 0.074 0.278 1.036 54.815 61.168
gp_m 0.700 0.700 0.010 0.010 0.683 0.716 1.002 2132.398 1884.079
Code
summ2 <- fit_rec$summary()

summ2[1:4,] %>% kable(digits = 3, caption = "Recovery")
Recovery
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -651.882 -651.798 10.587 10.281 -668.996 -634.343 1.003 638.313 1121.272
gp_theta 0.327 0.327 0.103 0.102 0.159 0.495 1.000 4799.380 3183.915
gp_sigma 0.754 0.750 0.114 0.115 0.571 0.947 1.001 996.214 1169.174
gp_m 0.086 0.074 0.063 0.063 0.007 0.203 1.000 3872.227 2602.156
Code
merged_df <- merge_prior_post(fit_gen,fit_rec,gp_theta, gp_sigma,gp_m)

merged_df %>% 
  # filter(model == "prior") %>% 
 ggplot(aes(y = .variable, x = .value)) + 
  stat_dotsinterval(aes(group = model, 
                        slab_colour = model, 
                        slab_fill = model,
                        interval_colour = model,
                        point_colour = model)) +
  facet_wrap(~.variable, scales = "free_x") +
  ggtitle("PPC of group level parameters")

Code
#ppc
fit_rec %>% spread_draws(subj_theta[n]) %>% ggplot(aes(x = subj_theta)) + stat_dots(quantiles = 100, fill = "red", alpha = 0.5) + stat_dots(data = as.data.frame(theta), aes(x = theta), quantiles = 100, fill = "blue", alpha = 0.5, inherit.aes = TRUE) + ggtitle("PPC of subject level theta parameters, Red = generated/estimated, Blue = real/ground truth")

Code
fit_rec %>% spread_draws(y_pred[n,t]) %>% sample_draws(ndraws = 1, seed = 789) %>% ggplot(aes(x = y_pred)) + stat_dots(quantiles = 100, fill = "red", alpha = 0.5) + stat_dots(data = as.data.frame(as.vector(y)), aes(x = `as.vector(y)`+0.1), quantiles = 100, fill = "blue", alpha = 0.5, inherit.aes = TRUE) + ggtitle("PPC Y, Red = generated/estimated, Blue = real/ground truth") + ggtitle("PPC of Bernoulli data, Red = generated/estimated, Blue = real/ground truth (shifted right to see easier)")

This only makes the posterior estimates worse!

Code
mod$code()
 [1] "data {"                                                                 
 [2] "  int<lower=0> N;"                                                      
 [3] "  int<lower=0> T;"                                                      
 [4] "  array[N,T] int y;"                                                    
 [5] ""                                                                       
 [6] "  real gp_theta_mu_pr;"                                                 
 [7] "  real gp_theta_sigma_pr;"                                              
 [8] "  "                                                                     
 [9] "  real gp_sigma_mu_pr;"                                                 
[10] "  real gp_sigma_sigma_pr;"                                              
[11] ""                                                                       
[12] "  real gp_m_mu_pr;"                                                     
[13] "  real gp_m_sigma_pr;"                                                  
[14] "  "                                                                     
[15] "  vector[N] psi_y;"                                                     
[16] "}"                                                                      
[17] ""                                                                       
[18] "parameters {"                                                           
[19] "  real gp_theta;"                                                       
[20] "  real<lower=0> gp_sigma;"                                              
[21] "  real<lower=0,upper=1> gp_m;"                                          
[22] ""                                                                       
[23] "  vector[N] subj_theta;"                                                
[24] ""                                                                       
[25] "}"                                                                      
[26] ""                                                                       
[27] "model{"                                                                 
[28] "    gp_theta ~ normal(gp_theta_mu_pr,gp_theta_sigma_pr);"               
[29] "    gp_sigma ~ normal(gp_sigma_mu_pr,gp_sigma_sigma_pr);"               
[30] "    gp_m ~ normal(gp_m_mu_pr,gp_m_sigma_pr);"                           
[31] "    "                                                                   
[32] "        for (n in 1:N){"                                                
[33] "          subj_theta[n] ~ normal(gp_theta + gp_m * psi_y[n], gp_sigma);"
[34] "        for (t in 1:T){"                                                
[35] "        y[n,t] ~ bernoulli_logit(subj_theta[n]);"                       
[36] "}"                                                                      
[37] "}"                                                                      
[38] "}"                                                                      
[39] ""                                                                       
[40] "generated quantities{"                                                  
[41] "  array[N,T] int y_pred;"                                               
[42] "  array[N,T] real log_lik;"                                             
[43] "  "                                                                     
[44] "          "                                                             
[45] "      for (n in 1:N){"                                                  
[46] ""                                                                       
[47] "      for (t in 1:T){"                                                  
[48] "        y_pred[n,t] = bernoulli_logit_rng(subj_theta[n]);"              
[49] "        "                                                               
[50] "        log_lik[n,t] = bernoulli_logit_lpmf(y[n,t] | subj_theta[n]);"   
[51] "       "                                                                
[52] "}"                                                                      
[53] "}"                                                                      
[54] "}"